      FUNCTION SIN_SUN(imonth,iday,ihour,iminu,iseco,xlatit)
!-----Computes sine of sun height angle (cosine of zenith angle)---------   
      implicit none
      
      real(8), intent(in) :: xlatit 
      integer(4), intent(in) :: imonth,iday,ihour,iminu,iseco

      real(8) :: SIN_SUN

      real(8) :: month,day,hour,minute,second,ndaym(12), &
      & days,delta,pi,tha,phi1
      integer(4) :: i
      data ndaym /31.,28.,31.,30.,31.,30.,31., &
      &           31.,30.,31.,30.,31./
      

      pi = 4.*atan(1.)
      phi1 = xlatit*pi/180.
      month = float(imonth)
      day   = float(iday)
      hour  = float(ihour)
      minute = float(iminu)
      second = float(iseco)
           
      if (month/=1.) then
        days=0.
        do i = 1, int(month)-1
          days = days + ndaym(i)
        enddo
        days = days + day
      else
        days = day
      endif
      days = days + int(hour/24.)

      delta = 23.5*pi/180.*dCOS(2*pi*(days-173.)/365.)
      tha   = 2*pi*(hour+minute/60.+second/3600.-12.)/24.
      SIN_SUN=dSIN(phi1)*dSIN(delta)+dCOS(phi1)*dCOS(delta)*dCOS(tha)

      END FUNCTION SIN_SUN
